﻿using System;

namespace LxCFD
{
    public class NSEquationBase : EquationBase
    {
        public Value1Base<double> Conv;//流动项前面的系数，动量方程为1，对于能量方程为比热容Cp
        public Value1Base<double> Gamma;//Tau由Continuum设置，但在各Zone中值有差异，因此另存一份，其他如源项则不再单独存储
        protected double[] MassImbalance;
        protected double[] MassFlux;
        protected float[] VelocityDirectionCorrection;//修正方向,见笔记
        protected ICalgCfPhir CalgCfPhir;
        #region 内部被调用的计算方法
        protected void set_x_flow(Value1Base<double> baseValue, IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            int ncell = 0;
            if (i1 < endp1) ncell = Faces[i1].NSideCell;
            for (; i1 < endp1; i1++, ncell++)
            {
                if (MassFlux[i1] > 0)//由明确的流出时，确定为流出边界
                {
                    CalBoundValue_grad0(ref Faces[i1]);
                }
                else
                {
                    x[ncell] = baseValue.GetCellValue(ncell);
                }
            }
        }
        protected double calUpwindValue_faceFluxPositive(ref FvFace face, double phiO, double phiN, double velDirCor)
        {
            double upwindvalue = phiN - 2 * Vector3d.DotProduct(Grad[face.OSideCell], face.OToN);
            return upwindvalue * velDirCor + phiO * (1 - velDirCor);
        }
        protected double calUpwindValue_faceFluxNegative(ref FvFace face, double phiO, double phiN, double velDirCor)
        {
            double upwindvalue = phiO + 2 * Vector3d.DotProduct(Grad[face.NSideCell], face.OToN);
            return upwindvalue * velDirCor + phiN * (1 - velDirCor);
        }

        protected void set_Ab_Upwind_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                double taugdiff = gamma * face.gDiff;
                double massflux = MassFlux[i1] * Conv.GetFaceValue(i1);
                b_FaceCommon[i1] = gamma * face.CalNonOrthogonalityCorrection_inner(Grad, x);
                if (massflux >= 0)
                {
                    ALUO[i1] = -taugdiff;
                    ALUN[i1] = -(massflux + taugdiff);
                }
                else
                {
                    ALUO[i1] = massflux - taugdiff;
                    ALUN[i1] = -taugdiff;
                }
            }
        }
        protected void set_Ab_SOU_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                double conv = Conv.GetFaceValue(i1);
                double gamadiff = gamma * face.gDiff;
                double massflux = MassFlux[i1] * conv;
                double facecor = gamma * face.CalNonOrthogonalityCorrection_inner(Grad, x);
                double veldircor = Math.Exp(-VelocityDirectionCorrection[i1] * conv / gamma);
                if (massflux >= 0)
                {
                    double phiOO = calUpwindValue_faceFluxPositive(ref face, x[face.OSideCell], x[face.NSideCell], veldircor);
                    double massgo = massflux * face.gOwner;
                    b_FaceCommon[i1] = massgo * phiOO + facecor;
                    ALUO[i1] = -gamadiff;
                    ALUN[i1] = -(massflux + massgo + gamadiff);
                }
                else
                {
                    double phiNN = calUpwindValue_faceFluxNegative(ref face, x[face.OSideCell], x[face.NSideCell], veldircor);
                    double massgn = massflux * (1 - face.gOwner);
                    b_FaceCommon[i1] = massgn * phiNN + facecor;
                    ALUO[i1] = massflux + massgn - gamadiff;
                    ALUN[i1] = -gamadiff;
                }
            }
        }
        protected void set_Ab_OSHER_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                double conv = Conv.GetFaceValue(i1);
                double gamadiff = gamma * face.gDiff;
                double massflux = MassFlux[i1] * conv;
                double phiO = x[face.OSideCell];
                double phiN = x[face.NSideCell];
                double facecor = gamma * face.CalNonOrthogonalityCorrection_inner(Grad, x);
                double veldircor = Math.Exp(-VelocityDirectionCorrection[i1] * conv / gamma);
                if (massflux >= 0)
                {
                    double phiOO = calUpwindValue_faceFluxPositive(ref face, phiO, phiN, veldircor);
                    double N_O = phiN - phiO;
                    if (N_O != 0)
                    {
                        double r = (phiO - phiOO) / N_O;
                        if (r < 0)
                        {
                        }
                        else if (face.gOwner * r < 1)
                        {
                            double massgo = massflux * face.gOwner;
                            b_FaceCommon[i1] = massgo * phiOO + facecor;
                            ALUO[i1] = -gamadiff;
                            ALUN[i1] = -(massflux + massgo + gamadiff);
                            continue;
                        }
                        else
                        {
                            b_FaceCommon[i1] = -massflux * N_O + facecor;
                            ALUO[i1] = -gamadiff;
                            ALUN[i1] = -(massflux + gamadiff);
                            continue;
                        }
                    }
                    b_FaceCommon[i1] = facecor;
                    ALUO[i1] = -gamadiff;
                    ALUN[i1] = -(massflux + gamadiff);
                    continue;
                }
                else
                {
                    double phiNN = calUpwindValue_faceFluxNegative(ref face, phiO, phiN, veldircor);
                    double O_N = phiO - phiN;
                    if (O_N != 0)
                    {
                        double gc = 1 - face.gOwner;
                        double r = (phiN - phiNN) / O_N;
                        if (r < 0)
                        {
                        }
                        else if (gc * r < 1)
                        {
                            double massgn = massflux * gc;
                            b_FaceCommon[i1] = massgn * phiNN + facecor;
                            ALUO[i1] = massflux + massgn - gamadiff;
                            ALUN[i1] = -gamadiff;
                            continue;
                        }
                        else
                        {
                            b_FaceCommon[i1] = -massflux * O_N + facecor;
                            ALUO[i1] = massflux - gamadiff;
                            ALUN[i1] = -gamadiff;
                            continue;
                        }
                    }
                    b_FaceCommon[i1] = facecor;
                    ALUO[i1] = massflux - gamadiff;
                    ALUN[i1] = -gamadiff;
                    continue;
                }
            }
        }
        protected void set_Ab_MUSCL_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                double conv = Conv.GetFaceValue(i1);
                double taugdiff = gamma * face.gDiff;
                double massflux = MassFlux[i1] * conv;
                double phiO = x[face.OSideCell];
                double phiN = x[face.NSideCell];
                double facecor = gamma * face.CalNonOrthogonalityCorrection_inner(Grad, x);
                double veldircor = Math.Exp(-VelocityDirectionCorrection[i1] * conv / gamma);
                if (massflux >= 0)
                {
                    double phiOO = calUpwindValue_faceFluxPositive(ref face, phiO, phiN, veldircor);
                    double N_O = phiN - phiO;
                    if (N_O != 0)
                    {
                        double r = (phiO - phiOO) / N_O;
                        if (r < LxTool.OneThird)
                        {
                            if (r < 0)
                            {
                            }
                            else
                            {
                                double mass2go = massflux * 2 * face.gOwner;
                                b_FaceCommon[i1] = mass2go * phiOO + facecor;
                                ALUO[i1] = -taugdiff;
                                ALUN[i1] = -(massflux + mass2go + taugdiff);
                                continue;
                            }
                        }
                        else
                        {
                            if (r * face.gOwner < 2 - face.gOwner)
                            {
                                double mass05go = massflux * 0.5 * face.gOwner;
                                b_FaceCommon[i1] = mass05go * phiOO + facecor;
                                ALUO[i1] = mass05go - taugdiff;
                                ALUN[i1] = -(massflux + taugdiff);
                                continue;
                            }
                            else
                            {
                                b_FaceCommon[i1] = -massflux * N_O + facecor;
                                ALUO[i1] = -taugdiff;
                                ALUN[i1] = -(massflux + taugdiff);
                                continue;
                            }
                        }
                    }
                    b_FaceCommon[i1] = facecor;
                    ALUO[i1] = -taugdiff;
                    ALUN[i1] = -(massflux + taugdiff);
                    continue;
                }
                else
                {
                    double phiNN = calUpwindValue_faceFluxNegative(ref face, phiO, phiN, veldircor);
                    double O_N = phiO - phiN;
                    if (O_N != 0)
                    {
                        double gc = 1 - face.gOwner;
                        double r = (phiN - phiNN) / O_N;
                        if (r < LxTool.OneThird)
                        {
                            if (r < 0)
                            {
                            }
                            else
                            {
                                double mass2gn = massflux * 2 * gc;
                                b_FaceCommon[i1] = mass2gn * phiNN + facecor;
                                ALUO[i1] = massflux + mass2gn - taugdiff;
                                ALUN[i1] = -taugdiff;
                                continue;
                            }
                        }
                        else
                        {
                            if (r * gc < 2 - gc)
                            {
                                double mass05gn = massflux * 0.5 * gc;
                                b_FaceCommon[i1] = mass05gn * phiNN + facecor;
                                ALUO[i1] = massflux - taugdiff;
                                ALUN[i1] = -(mass05gn + taugdiff);
                                continue;
                            }
                            else
                            {
                                b_FaceCommon[i1] = -massflux * O_N + facecor;
                                ALUO[i1] = massflux - taugdiff;
                                ALUN[i1] = -taugdiff;
                                continue;
                            }
                        }
                    }
                    b_FaceCommon[i1] = facecor;
                    ALUO[i1] = massflux - taugdiff;
                    ALUN[i1] = -taugdiff;
                    continue;
                }
            }
        }
        protected void set_Ab_TVD_DC_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                double conv = Conv.GetFaceValue(i1);
                double taugdiff = gamma * face.gDiff;
                double massflux = MassFlux[i1] * conv;
                double phiO = x[face.OSideCell];
                double phiN = x[face.NSideCell];
                double facecor = gamma * face.CalNonOrthogonalityCorrection_inner(Grad, x);
                double veldircor = Math.Exp(-VelocityDirectionCorrection[i1] * conv / gamma);
                if (massflux >= 0)
                {
                    double N_O = phiN - phiO;
                    double r = 1;
                    if (N_O != 0)
                    {
                        r = 2 * Vector3d.DotProduct(Grad[face.OSideCell], face.OToN) / N_O - 1;
                        r *= veldircor;
                        r = CalgCfPhir.cal(r, face.gOwner);
                    }
                    b_FaceCommon[i1] = -massflux * r * N_O + facecor;
                    ALUO[i1] = -taugdiff;
                    ALUN[i1] = -(massflux + taugdiff);
                }
                else
                {
                    double O_N = phiO - phiN;
                    double r = 1;
                    if (O_N != 0)
                    {
                        r = -2 * Vector3d.DotProduct(Grad[face.NSideCell], face.OToN) / O_N - 1;
                        r *= veldircor;
                        r = CalgCfPhir.cal(r, 1 - face.gOwner);
                    }
                    b_FaceCommon[i1] = -massflux * r * O_N + facecor;
                    ALUO[i1] = massflux - taugdiff;
                    ALUN[i1] = -taugdiff;
                }
            }
        }
        protected void set_Ab_boundOwner_DiffPart(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                ALUO[i1] = ALUN[i1] = -gamma * face.gDiff;//A_LU_N赋值是为了给face.ownercell用的，用来计算A_Diag
                b_FaceCommon[i1] = gamma * face.CalNonOrthogonalityCorrection_bound(Grad, x);
            }
        }
        protected void set_A_boundOwner_0(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ALUO[i1] = ALUN[i1] = 0;//A_LU_N赋值是为了给face.ownercell用的，用来计算A_Diag
                b_FaceCommon[i1] = 0;//如果数组共用，可能有错误，需赋值
            }
        }
        protected void set_Ab_bound_flow(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double massflux = MassFlux[i1] * Conv.GetFaceValue(i1);
                if (massflux > 0)//由明确的流出时，确定为流出边界
                {
                    //向外流出，边界条件为：梯度为0，因此A_LU_O设置为0
                    ALUO[i1] = 0;
                    ALUN[i1] = -massflux;
                    b_FaceCommon[i1] = massflux * (x[face.OSideCell] - x[face.NSideCell]);//用downwind格式
                    //对于N，不计算，通过updatebound算出
                }
                else
                {
                    //流入，边界条件为：固定值，A_LU_O选择计算（如果仍然假设梯度为0会怎样？）
                    double gamma = Gamma.GetFaceValue(i1);
                    double taugdiff = gamma * face.gDiff;
                    b_FaceCommon[i1] = gamma * face.CalNonOrthogonalityCorrection_bound(Grad, x);
                    ALUO[i1] = massflux - taugdiff;
                    ALUN[i1] = -taugdiff;
                    //对于N，不计算，通过updatebound算出
                }
            }
        }
        protected void setAb_step2(Value1Base<double> source1, int taskIndex)
        {
            if (source1.IsEqual0())//sb is Value1Const sc && sc.number == 0
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    ref FvCell cell = ref Cells[i1];
                    double massimb = MassImbalance[i1] * Conv.GetCellValue(i1);
                    ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO) - massimb;//有个规律，ADiag恰好是O面加上ALUN，N面加上ALUO
                    b[i1] = cell.CalIfOwnerAddElseMunus(b_FaceCommon);
                }
            }
            else
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    ref FvCell cell = ref Cells[i1];
                    double massimb = MassImbalance[i1] * Conv.GetCellValue(i1);
                    ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO) - massimb;//有个规律，ADiag恰好是O面加上ALUN，N面加上ALUO
                    b[i1] = cell.CalIfOwnerAddElseMunus(b_FaceCommon) + source1.GetCellValue(i1) * cell.Volumn;
                }
            }
        }

        protected void setAb_step2(Value2Base<double> source2, int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvCell cell = ref Cells[i1];
                double massimb = MassImbalance[i1] * Conv.GetCellValue(i1);
                source2.GetCellValue(out double S_Ad, out double S_b, i1);
                ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO) - massimb - S_Ad * cell.Volumn;
                b[i1] = cell.CalIfOwnerAddElseMunus(b_FaceCommon) + S_b * cell.Volumn;
            }
        }
        #endregion
        protected interface ICalgCfPhir
        {
            double cal(double r, double gCf);//计算结果位置为phiC+（结果）*(PhiD-PhiC)
        }
        protected class MUSCL_DC : ICalgCfPhir
        {
            public double cal(double r, double gCf)
            {
                return Math.Max(0, Math.Min(Math.Min(2 * gCf * r, (r + 1) * gCf / 2), 1));
                // return Math.Max(0, Math.Min(Math.Min(4 * gCf * r, (r + 1) * gCf), 2));
            }
        }
        public class OSHER_DC : ICalgCfPhir
        {
            public double cal(double r, double gCf)
            {
                return Math.Max(0, Math.Min(gCf * r, 1));
            }
        }
    }
}
//public override void SetAb_Step2(int taskIndex)
//{
//    if (Source is Value1Base sb)
//    {
//        if (sb is Value1Const sc && sc.number == 0)
//        {
//            InnerIndex.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//            for (; i1 < endp1; i1++)
//            {
//                ref Mesh.Cell cell = ref cells[i1];
//                double massimb = massImbalance[i1] * Conv.getCellValue(i1);
//                if (massimb >= 0)//流出的多流入的少，修正时假设引入mass，标量值按上次迭代
//                {
//                    ADiag[i1] = -cell.calAllAdd(ALUN, ALUO);//有个规律，ADiag恰好是O面加上ALUN，N面加上ALUO
//                    b[i1] = cell.calIfOwnerAddElseMunus(b_faceCommon) + massimb * x[i1];
//                }
//                else//流出的少，流入的多，修正时假设流出部分mass，流出的参数按本网格值，
//                {
//                    ADiag[i1] = -cell.calAllAdd(ALUN, ALUO) - massimb;//有个规律，ADiag恰好是O面加上ALUN，N面加上ALUO
//                    b[i1] = cell.calIfOwnerAddElseMunus(b_faceCommon);
//                }
//            }
//        }
//        else
//        {
//            InnerIndex.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//            for (; i1 < endp1; i1++)
//            {
//                ref Mesh.Cell cell = ref cells[i1];
//                double massimb = massImbalance[i1] * Conv.getCellValue(i1);
//                if (massimb >= 0)//流出的多流入的少，修正时假设引入mass，标量值按上次迭代
//                {
//                    ADiag[i1] = -cell.calAllAdd(ALUN, ALUO);//有个规律，ADiag恰好是O面加上ALUN，N面加上ALUO
//                    b[i1] = cell.calIfOwnerAddElseMunus(b_faceCommon) + massimb * x[i1] + sb.getCellValue(i1) * cell.Volumn;
//                }
//                else//流出的少，流入的多，修正时假设流出部分mass，流出的参数按本网格值，
//                {
//                    ADiag[i1] = -cell.calAllAdd(ALUN, ALUO) - massimb;//有个规律，ADiag恰好是O面加上ALUN，N面加上ALUO
//                    b[i1] = cell.calIfOwnerAddElseMunus(b_faceCommon) + sb.getCellValue(i1) * cell.Volumn;
//                }
//            }
//        }
//    }
//    else if (Source is Value2Base sb2)
//    {
//        InnerIndex.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//        for (; i1 < endp1; i1++)
//        {
//            ref Mesh.Cell cell = ref cells[i1];
//            double massimb = massImbalance[i1] * Conv.getCellValue(i1);
//            sb2.getCellValue(out double S_Ad, out double S_b, i1);
//            if (massimb >= 0)//流出的多流入的少，修正时假设引入mass，标量值按上次迭代
//            {
//                ADiag[i1] = -cell.calAllAdd(ALUN, ALUO) - S_Ad * cell.Volumn;
//                b[i1] = cell.calIfOwnerAddElseMunus(b_faceCommon) + massimb * x[i1] + S_b * cell.Volumn;
//            }
//            else//流出的少，流入的多，修正时假设流出部分mass，流出的参数按本网格值，
//            {
//                ADiag[i1] = -cell.calAllAdd(ALUN, ALUO) - massimb - S_Ad * cell.Volumn;
//                b[i1] = cell.calIfOwnerAddElseMunus(b_faceCommon) + S_b * cell.Volumn;
//            }
//        }
//    }
//    else
//    {
//        LxConsole.AddLine(LogTypes.错误, "未知的源项类型");
//    }
//}